Semivariogramas univariados
vgmno2 <- vgm(psill = mat1[1, 1],
model = "Sph",
range = 6096,
add.to = vgm(psill = mat2[1, 1],
model = "Hol",
range = 2294))
vgmo3 <- vgm(psill = mat1[2, 2],
model = "Sph",
range = 6096,
add.to = vgm(psill = mat2[2, 2],
model = "Hol",
range = 2294))
vgmnox <- vgm(psill = mat1[3, 3],
model = "Sph",
range = 6096,
add.to = vgm(psill = mat2[3, 3],
model = "Hol",
range = 2294))
Semivarogramas cruzados (Bivariados)
vgmno2_o3 <- vgm(psill = mat1[1, 2], model = "Sph",
range = 6096,
add.to = vgm(psill = mat2[1, 2],
model = "Hol",
range = 2294))
vgmno2_nox <- vgm(psill = mat1[1, 3],
model = "Sph",
range = 6096,
add.to = vgm(psill = mat2[1, 3],
model = "Hol",
range = 2294))
vgmno3_nox <- vgm(psill = mat1[2, 3],
model = "Sph",
range = 6096,
add.to = vgm(psill = mat2[2, 3],
model = "Hol",
range = 2294))
gstat
remove_na <- function(frame, vari_) {
# Remove na from sp object
datos1 <- frame
bool <- !is.na(datos1@data[vari_])
datos1@data <- datos1@data[bool, ]
datos1@coords <- datos1@coords[bool, ]
return(datos1)
}
coordinates(datos) <- ~ X + Y
g_st <- gstat(NULL,
id = "NO2",
formula = NO2 ~ X + Y,
model = vgmno2,
data = remove_na(datos, "NO2"))
g_st <- gstat(g_st,
id = "O3",
formula = O3 ~ Y,
model = vgmo3,
data = remove_na(datos, "O3"))
g_st <- gstat(g_st,
id = "NOX",
formula = NOX ~ Y,
model = vgmnox,
data = remove_na(datos, "NOX"))
#Cruzados
g_st <- gstat(g_st,
id = c("NO2", "O3"),
model = vgmno2_o3)
g_st <- gstat(g_st,
id = c("NO2", "NOX"),
model = vgmno2_nox)
g_st <- gstat(g_st,
id = c("O3", "NOX"),
model = vgmno3_nox)
pander::pander(do.call(rbind, g_st$model)[, 1:3])
| NO2.1 |
Hol |
13.02 |
2294 |
| NO2.2 |
Sph |
30 |
6096 |
| O3.1 |
Hol |
46.4 |
2294 |
| O3.2 |
Sph |
50 |
6096 |
| NOX.1 |
Hol |
26.96 |
2294 |
| NOX.2 |
Sph |
35 |
6096 |
| NO2.O3.1 |
Hol |
24.54 |
2294 |
| NO2.O3.2 |
Sph |
30 |
6096 |
| NO2.NOX.1 |
Hol |
18.73 |
2294 |
| NO2.NOX.2 |
Sph |
30 |
6096 |
| O3.NOX.1 |
Hol |
35.36 |
2294 |
| O3.NOX.2 |
Sph |
30 |
6096 |
Mapas de predicción de O3 con las covariables espaciales NO2 y NOX
prediction_plot <- function(g_object, variable, map_path) {
map <- readOGR(map_path)
new <- sp::spsample(map, n = 100000, type = "regular")
coordinates(new) ~ x1 + x2
colnames(new@coords) <- c("X", "Y")
predic <- predict(g_object, newdata = new)
prediction <- data.frame(predic)
pred <- paste(variable, ".pred", sep = "")
plot <- ggplot(prediction, aes_string("X", "Y", fill = pred)) +
geom_tile() +
scale_fill_viridis_c() +
theme_void()
return(plot)
}
variance_plot <- function(g_object, variable, map_path) {
map <- readOGR(map_path)
new <- sp::spsample(map, n = 10000, type = "regular")
coordinates(new) ~ x1 + x2
colnames(new@coords) <- c("X", "Y")
predic <- predict(g_object, newdata = new)
prediction <- data.frame(predic)
var <- paste(variable, ".var", sep = "")
plot <- ggplot(prediction, aes_string("X", "Y", fill = var)) +
geom_tile() +
scale_fill_viridis_c(option = "inferno",
direction = -1) +
theme_void()
return(plot)
}
cv_plot <- function(g_object, variable, map_path) {
map <- readOGR(map_path)
new <- sp::spsample(map, n = 10000, type = "regular")
coordinates(new) ~ x1 + x2
colnames(new@coords) <- c("X", "Y")
predic <- predict(g_object, newdata = new)
prediction <- data.frame(predic)
pred <- paste(variable, ".pred", sep = "")
var <- paste(variable, ".var", sep = "")
aux <- abs(sqrt(prediction[var]) / abs(prediction[pred]))
aux[aux > 1] <- 1
prediction["cv"] <- aux
plot <- ggplot(prediction, aes_string("X", "Y", fill = "cv")) +
geom_tile() +
scale_fill_viridis_c(option = "magma",
direction = -1) +
theme_void()
return(plot)
}
pl1 <- prediction_plot(g_st, "O3",
"2_COK_G_stat/SP/mpiosutm.shp")
## OGR data source with driver: ESRI Shapefile
## Source: "/home/jncc/Documents/Notebooks/Monitorias/Bookdown/Cuadernos/2_COK_G_stat/SP/mpiosutm.shp", layer: "mpiosutm"
## with 54 features
## It has 7 fields
## Linear Model of Coregionalization found. Good.
## [using universal cokriging]
pl2 <- variance_plot(g_st, "O3",
"2_COK_G_stat/SP/mpiosutm.shp")
## OGR data source with driver: ESRI Shapefile
## Source: "/home/jncc/Documents/Notebooks/Monitorias/Bookdown/Cuadernos/2_COK_G_stat/SP/mpiosutm.shp", layer: "mpiosutm"
## with 54 features
## It has 7 fields
## Linear Model of Coregionalization found. Good.
## [using universal cokriging]
pl3 <- cv_plot(g_st, "O3",
"2_COK_G_stat/SP/mpiosutm.shp")
## OGR data source with driver: ESRI Shapefile
## Source: "/home/jncc/Documents/Notebooks/Monitorias/Bookdown/Cuadernos/2_COK_G_stat/SP/mpiosutm.shp", layer: "mpiosutm"
## with 54 features
## It has 7 fields
## Linear Model of Coregionalization found. Good.
## [using universal cokriging]